Accurate extraction of the News 
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We propose a new scheme for extracting gravitational radiation from a characteristic numerical 
simulation of a spacetime. This method is similar in conception to our earlier work but analytical 
and numerical implementation is different. The scheme is based on direct transformation to the 
Bondi coordinates and the gravitational waves are extracted by calculating the Bondi news function 
in Bondi coordinates. The entire calculation is done in a way which will make the implementation 
easy when we use uniform Bondi angular grid at T"*" . Using uniform Bondi grid for news calculation 
has added advantage that we have to solve only ordinary differential equations instead of partial 
differential equation. For the test problems this new scheme allows us to extract gravitational 
radiation much more accurately than the previous schemes. 

PACS numbers: 04.25.Dm, 04.30.-w, 95.85.Sz 



I. INTRODUCTION 



One of the main aims of the numerical relativity sim- 
ulations of spacetimes is to calculate the gravitational 
radiation emitted by various kinds of sources, particu- 
larly in the strong field regime. The predicted waveforms 
for different kinds of astrophysical sources will be useful 
for detection of gravitational waves and parameter esti- 
mation of sources (by matched filtering technique) using 
many working and upcoming gravitational wave detec- 
tors hke, LIGO, VIRGO, TAMA, GEO and LISA. For 
this we need the simulations to be very accurate and we 
also need to extract the gravitational waves faithfully and 
accurately from the simulations. Both these tasks are 
very challenging. Binary black holes and neutron stars 
are the main likely sources for early detection. In the re- 
cent times there is considerable progress in evolving the 
binary black hole Q, S H 0, S |SB II 13 and binary 
neutron star [10| | spacetimes, but a lot of improvements 
are needed for the predicted waveforms to be useful for 
the future observations. 

In the most common approach to the numerical relativ- 
ity (Cauchy/ADM approach) the spacetime is foliated by 
a sequence of spacelike hypersurfaces parameterized by 
time. In such a case we cannot extract the gravitational 
waves exactly and some approximate techniques have to 
be developed. These techniques [H El El El El [11 
use the data at the outer boundary of the simulations to 
predict the gravitational wave content of spacetimes at 
future null infinity {l^'^). In recent times new promising 
ways, which are also useful for extracting gravitational 
waves from Cauchy simulations of spacetimes have been 
proposed [H El 111 HI IH IH . But in this paper we 
use the characteristic (null cone) formulation using the 
Bondi-Sachs coordinates. The spacetime is foliated by 
series of null hypersurfaces parameterized by outgoing 
(/ingoing) null geodesies. Using compactification tech- 



1shrir@hri.res.ml 



niques, X"*" is included in the computational domain. In 
characteristic numerical simulations one can in principle 
extract the gravitational radiation exa ctly at null infinity 
(X+) using the Bondi news function |3 HI IH- Even 
then the accurate extraction of the news seems quite 
tricky as many numerical complications make the task 
very hard. The news function takes a very simple form 
in the Bondi coordinates, but the numerical simulations 
use more general Bondi-Sachs coordinates. The Bondi 
coordinates correspond to inertial observers at but 
in a numerical code we want to provide boundary con- 
ditions at inner boundary (e.g. black hole horizon) and 
then solve the equations as we go out to X'^ , so we use 
the more general Bondi-Sachs coordinates. Earlier news 
was computed directly in Bondi-Sachs coordinates 12^ 
and then a transformation was done to the Bondi coor- 
dinates (as detectors are more or less in the reference 
frame of inertial observer at I^). The news expression 
in Bondi-Sachs coordinates is quite complicated. 

We had proposed a scheme [23 where we calculated the 
news by going directly to the Bondi coordinates. But fur- 
ther tests showed that improvement in the news calcula- 
tions were needed. In certain cases numerically interme- 
diate constraint, Jo = (it basically says, the tangential 
part of Bondi metric at X+ should be zero, the notation is 
introduced in Sec. was not getting satisfied, though as 
such the final news calculations gave satisfactory result. 

The problem was traced to a particular (i.e. dX) term 
in equation (20) for Jq in [2l|. This term is basically 
related to an angular derivative of the Bondi angular 
coordinate with respect to the Bondi-Sachs angular co- 
ordinate. Many techniques, like doing entire numerical 
calculations using spin weighted spherical harmonic co- 
efficients etc. were tried to improve the accuracy of the 
result while sticking to the original analytical procedure 
in '2^, but it was hard to achieve the desired results. 

The main reason for this problem seems to be the in- 
teraction of various numerical errors creating unexpected 
enhancement of error while calculating angular deriva- 
tives. Using the scheme in 28], these angular deriva- 
tives had to be calculated in a complicated way, using 
Jacobian of coordinate transformation. This was needed 
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as all the angular derivatives were with respect to the 
Bondi-Sachs coordinates while for calculating the news 
we needed to use uniform Bondi angular grid. This choice 
had to made so that one could interpolate things more 
easily and also as then one needed to solve only ordinary 
differential equations for doing the coordinate transfor- 
mation instead of partial differential equations. 

In this paper we present a scheme which suites these 
numerical requirements better and gives the desired accu- 
racy. The scheme is similar in conception to the scheme 
in |28j | , but the analytical and numerical implementation 
is different. The news is calculated by directly going 
to the Bondi coordinates. At the initial slice Jq has to 
be zero (notation is introduced in section Sec. and 
the Bondi-Sachs and Bondi coordinates have to match 
at X+. For the news calculation we stick to the Bondi 
coordinates and uniform Bondi grid at X+ . The transfor- 
mation is done analytically (partly using Mathematica) 
so that the numerical implementation will have good ac- 
curacy. We have implemented our scheme and tested it 
with analytical solutions. The calculations and results 
are presented in this paper. This new scheme shows very 
good accuracy for the news extraction and also the con- 
straints Jo = is satisfied with desired accuracy. 

The paper starts with a summary of relevant results 
and notation for the characteristic formulation of numer- 
ical relativity (Sec. |nj. The coordinate transformation 
and algebraic calculations and expressions are given in 
Sec. mil and then the procedure for computing the news, 
at both analytic and computational levels, is described 
in Sec. IIVI The computational tests and results are pre- 
sented in Sec. followed by discussions and conclusions 
in Sec. ED 



II. NOTATION 

Here we briefly review the notation and formalism 
for the characteristic numerical relativity [23, (see 
also 0, Isol I31I Is^ ) . The formalism is based on a fam- 
ily of outing null hypersurfaces using the Bondi-Sachs 
coordinates. The hypersurfaces are labeled by u, 
{A = 2,3) label the null rays and r is a surface area 
coordinate. The metric in these Bondi-Sachs coordi- 
nates I23, Isal is written as, 



e^P{\ + —)-r^hABU''U 
r 



-2e^l^dudr ~ 2r'^hARU^ dudx 



r'hABdx'^dx'', (1) 

where h^^hBC ~ det{hAB) — det{qAB), with 

QAB a unit sphere metric. We use stereographic coor- 
dinates x^ = {q,p) for which the unit sphere metric is 



qABdx'^dx^ = -^{dq^ + dp^), 



where 



P^l + q^+p^ 



(2) 
(3) 



We introduce a complex dyad qA defined by 



(4) 



with i = For an arbitrary Bondi-Sachs metric, Hab 

can then be represented by its dyad component 



J 



hABq^q"" I'i. 



(5) 



hAB is uniquely determined by J, since the determinant 
condition [detiJiAB) — det{qAB)) implies that there are 
only two independent components of Hab and the re- 
maining dyad component 

K = hABq^f/2 (6) 

satisfies I — — J J. For spherically symmetric case J 
is identically zero. We introduce the spin-weighted field. 



U = U^qA. 



(7) 



We also introduce complex angular coordinate ( = q + ip 
as well as the (complex differential) eth operators 3 and 
§ (see [S^ for full details) which are given as. 



dA 
9C 



sAC dA = P 



dA 



sAC, 



(8) 



where A is any spin weighted field with spin weight s. 

The news is calculated in a conformally compactified 
coordinates. Specifically, {u,r,x^) — > {u,£,x'^) where 
£= 1/r. In {u,£,x^) coordinates, the compactified met- 
ric is ds^ = f'ds^ where £ is a conformal factor with 
future null infinity given by ^ = 0. The compactified 
metric will be denoted by g"^ (a, = — 3) and the 
general compactified Bondi-Sachs metric is 



^11 



e'^^'Va, 9 



lA 



9'" = e 



-2/3 



t.AB 



= h 



AB 



9"^ = r - 0, 



(9) 



where Va = £'^{1 + £W). In addition to the general com- 
pactified Bondi-Sachs metric Eq. (|5J), we will refer to 
the compactified Bondi-Sachs metric satisfying the Bondi 
conditions, and such quantities will be denoted with a 
tilde {'), with compactified metric and coordinates g"^ 
and {u,£,x^), respectively. OnX+, i.e. £ ^0, 9°"^ satis- 
fies the Bondi conditions 

500 = 0, = 0, goi = 1, gAB^qAB, (10) 

where q'^^ is a unit sphere metric with respect to the 
Bondi angular coordinates i"^. 



III. COORDINATE TRANSFORMATION 

We define a coordinate transformation, near , be- 
tween (u, ^, x^) and {u, £, x^) as a Taylor series expansion 
in t. Up to second order in t we can write, 

u = uq + A^^i + C"P, £ = £/oj + C^?, 

x^ = x^ + A^i+CH^, (11) 
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where co, Xq, A-^, uq, are all functions of and u 
only. It will turn out that the C" are not needed for the 
news calculation, but they are included so that the Jaco- 
bian for coordinate transformation is manifestly correct 
to first order in L We also introduce complex quantities 

A = qAA^ and X = qax^- (12) 

The metric functions are expanded to first order in £ as, 



(13) 



The same first-order expansion in I and notation will be 
used for the metric quantities in the coordinate system 
satisfying the Bondi conditions (e.g., J = Jq +iJi). The 
general compactificd Bondi-Sachs metric, and the Bondi- 
Sachs metric satisfying the Bondi conditions, are related 

by 



{gap)o + K9af3)l = gap = 



, dxf^ dx" 



(14) 



95" 95/5' 

with the factor appearing because there is also an 
implicit change of compactification factor from £ to £. 

We use Mathematica to find g"^ to first order in £, and 
in doing so we have imposed the conditions, 

du 

du 

and 



(15) 



du 



2/3 



These conditions are equivalent of satisfying, 



{du + U^dB)i^ - 



(16) 



(17) 



and 



200 



(18) 

along the generators of 1+ (see Eqs. (39) and (40) in j27| 
for some details). This ensures Bondi conditions 500 = Oj 
goA = 0, goi = 1 are satisfied identically. 

We also evolve uj along the generators of 1+ by inte- 
grating [23, 

' BU). (19) 



du 
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Evolving Lj using this ordinary differential equation re- 
duces the errors due to angular derivatives as compared 
to the method used in j23|. Knowing the coordinate 
transformation and w, we can find the value of Jo (Jo 
in Bondi coordinates). 



Jo 



4P2 

-2PPdu{jQUo + KQUo)iSX - 
-2PPdudX{KoUo + JqUq) 
+2p'^KodX{dX + BxPX)), 



BxPX) 



(20) 



Where, Bx = The analytical value of Jo has to be 
zero as we are in Bondi coordinates and we can use the 
numerical value to estimate the accuracy of the coordi- 
nate transformation. 
Bondi-Sachs conditions 



r - 0, 

jives us from the lowest (O"*) order in £, 
Au — — ^oje^^duBu. 



(21) 



(22) 



First order in £ will give us but, it is not needed for 
news calculation. Bondi-Sachs condition, 



g 



0, 



(23) 



leads to. 



A 



^(((9X -I- BxPX)P - BuUoP)du + duQXP) 

. (24) 
from the lowest (O*'') order in £, and will give C from 1"* 
order, but it is not needed for our calculation. One can 
also get Au and A by solving for 511 — and goA — but 
we use 3°*^ = 0, g'^^ = as it simplifies the calculation. 

After knowing these quantities, we can calculate J^. 
The expression for J^ is quite long and is given in the 
equation ljAl|l in the appendix. 



IV. CALCULATING THE NEWS 

After knowing the coordinate transformations to re- 
quired order we can develop a systematic scheme for ex- 
tracting the news. First we describe the analytic aspects 
and then we discuss some details about the actual im- 
plementation in the program. It is assumed that the 
Bondi-Sachs metric is knows throughout the spacetime. 

A. Procedure 

1. We integrate equation (fTH|l to get x^ on a fixed 
Bondi grid. 

2. We evolve equation ((T^ for uj. 

3. Jo is evaluated form equation (|20l) and is used to 
monitor the numerical accuracy of the calculation. 

4. We find the time transformation using equation 

5. Au is calculated using equation H22(l and A from 
equation (|^ . 

6. We evaluate Jf using equation ljAl|l . 
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7. The news is extracted by differentiating Jj, 



IdJi 
2 du 



(25) 



B. Computational implementation 

Tlie news module lias been written to interface directly 
with the null gravity code in its current form |35l |. Thus 
we use a compactified radial coordinate x = r/{R + r), 
with R = I. There are points in the x direction in 
the range [0.5, 1] (corresponding to 1 < r < oo). We use 
stereographic angular coordinates with two patches with 
minimum patch overlap of 5 grid points. Both patches 
have Un grid points in each angular direction. 

In actual numerical implementation for the news calcu- 
lation, we use uniform {u, £,x^) grid instead of uniform 
{u, £,x^) to avoid having interpolation with respect to 
time. This means in place of equation (|15|l we solve its 
inverse, 



du 



(26) 



along the generators of X+ . This also means that in the 
derivative operators with Bondi coordinates, have to be 
replaced by, 



and 



d 



d_ 

du 

d 



-2/3o Q 

Lv du 



3-2/30 Q 



(27) 



dx^ " d 



l^- (28) 



We also introduce dU and dU as well as EX and BX as 
additional variables. This helps to reduce the numerical 
errors as typically news indirectly involves 3'''' or higher 
angular derivative of U and X (if the original coordi- 
nates are not already Bondi coordinates). Higher angular 
derivatives are known to give lot of problems numerically. 
Introducing these additional variables is one simple way 
of improving the numerical behavior. These variables 
should be introduced in the null evolution code as well. 
All integration and differentiation schemes are 2"'^ order, 
apart from for the first step, for which time integration 
is done by 1^* order scheme. 



V. COMPUTATIONAL TESTS AND RESULTS 

We implemented and tested our method for calculating 
the Bondi news function for Schwarzschild and linearized 
Robinson- Trautman in rotating as well as tumbling coor- 
dinates. The analytical solution on the discretized grid is 
given as the input to the news module. The news mod- 
ule as such sees it as input from an equivalent numerical 
evolution code. The same news module can be used with 
suitable null evolution code. 



Schwarzschild solution in rotating and tumbling 
coordinates 



We get the Schwarzschild solution in Bondi-Sachs ro- 
tating coordinates by transforming the standard angular 
coordinates by, (f) cf) + ku. Various metric coefficients 
at are written in this case as, 



Wn = Ws^ 0, /3w = /3s = 0, 
Jn — Js = 0, {Je)N = {■h)s = 0, 



Un = - 



M{1 + CnCn) 



Us = 



M(l+CsCs)' 



(29) 



where the subscripts N and S refer to north and south 
patch. 

When we have coordinate system which is rotating 
around an equatorial axis, we call it a tumbling coordi- 
nate system. In this case the relationship of the tumbling 
angular coordinates to the Bondi angular coordinates is 
given as. 



6 — arccos(cos 9 cos ku + sin 9 sin (f) sin ku) 



arctan 



cos ku sin 9 sin ( 



cos 6 sin KU 



sm u cos ( 



(30) 



(31) 



In these Bondi-Sachs coordinates only nonzero part of 
the metric at is U, which is given as, 



Un = iK 



1-C 



N 



(1 + C/vC/v)' 



Us — iK 



(l + CsCs)' 



(32) 



The convergence tests were performed for the following 
grid sizes for all the tests. 



(«) 


Aq = 


Ap 


(b) 


Aq = 


Ap 


ic) 


Aq^ 


Ap 


id) 


Aq^ 


Ap 



1/8, Ax = 
1/16, Aa; 



1/64, Am = 0.04 
= 1/128, Aw = 0.02 



For the highest resolution run we keep Ax = 1/64 due 
to memory limitations, but as such the results do not 
change with Ax. 

In both these cases the L2 norm of Jq shows 2nd (1.99) 
order convergence to zero as in figure Q and the conver- 
gence rate is essentially independent of u. 

The news function with our scheme remains identically 
zero. This is an ideal behavior and is a pleasant surprise. 
The news shows this behavior as A and Au are zero an- 
alytically and with our scheme and choice of variables 
they also remain zero numerically. From the expression 
for Jf one can see that it will be zero for our case. Due to 
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FIG. 1: L2 norm of || J0II2 for Schwarzschild in tumbling coor- 
dinates as function of u. The norm changes very slowly with u. 
Analytically Jo should be zero. The graph for Schwarzschild 
in rotating coordinates looks essentially identical, so we have 
not plotted it separately. 



this we can expect to extract out news accurately even if 
we happen to be in a rotating or tumbUng kind of back- 
ground coordinates and the amplitude of gravitational 
radiation is very small. 

The scheme in showed required convergence for 
Jo in rotating Schwarzschild case, but later it was re- 
alized it had problems while dealing with the tumbling 
Schwarzschild case. It was very hard to remove these nu- 
merical problems in that scheme. Also the error in the 
news in both these cases, though, convergent to 2"'^ or- 
der, they were not all that small, if one had to extract 
very weak gravitational wave signals in nontrivial coor- 
dinates. Both these issues lead us to develop the scheme 
presented above. 



convergence rate at u 


grids (a) & (b) 


(b) & (c) 


(c) & (d) 


4 


2.012 


2.002 


2.008 


8 


2.028 


2.021 


2.005 


12 


2.070 


2.013 


2.005 


16 


2.097 


2.012 


2.005 


20 


2.086 


2.011 


2.005 



TABLE I: Change of convergence rate with u for the L2 
norm of error in the news function for linearized Robinson- 
Trautman in rotating coordinates for different grid resolutions 
in equation 1331 



1. Rotating coordinates 

The metric components for the linearized Robertson 
Trautmann solution, when only Y20 term is present (i.e. 
only the Z = 2, m = term of spherical harmonics is 
present in the perturbation) can be written at 2"+ as [s^ , 



Wn = Ws 

Jn ~ Js 
{Ji)n 

Un 
Us 



-12A,ne-2-/^AJ- 



N 



(1 + CnCnY 



-12A2oe"2''/^'^Af 



C^5 



(1 + CsCs) 



2 ' 



N 



M(1 + CwCjv)' 
M{1 + CsCs) ' 



(34) 



where the subscripts N and S refer to north and south 
patch, M gives the mass of the remnant black hole and 
A20 is the initial (at u = 0) amplitude of the perturbation. 

In this case Jq shows 2"'^ order convergence to zero. 
Basically the Jq part remains same as for the correspond- 
ing Schwarzschild case (Fig. QJ. The news also shows 2"*^ 
order convergence to the analytical value on the Bondi 
grid given by, 



B. Linearized Robinson- Trautman solution 

The Robinson- Trautman solution represents a dis- 
torted black hole emitting purely outgoing radiation. 
The radiation decays exponentially, and asymptotically 
the solution becomes Schwarzschild. In the linearized 
case when the amplitude of the perturbation (i.e. also 
the gravitational radiation) is small one can write it as 
addition of different spherical harmonics (Yim) of various 
amplitudes (see e.g. [23 )• 



N = -12A 



2oe 



(1 + CC)^ 



(35) 



on both the patches. 

For our test in figure we set M — I, A20 = 3 x 10"'' 
and K = 0.01. Even for very small or very hight value 
of A20, news shows 2"*^ order convergence though the 
gravitational wave amplitude goes down by more than 10 
order of magnitude as the time passes. The convergence 
rate are shown in table ((lT|) . This is much better behavior 
than we could get by the earlier news modules. 
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FIG. 2: L2 norm of the error in the news as a function of 
time u for various discretizations for the linearized Robinson- 
Trautman solution in rotating coordinates. 




FIG. 3: 1/2 norm of the error in the news as a function of 
time u for various discretizations for the hnearized Robinson- 
Trautman solution in tumbling coordinates. 



2. Tumbling coordinates 

Various metric components for the linearized Robert- 
son Trautmann solution, when only Y20 perturbation 
term is present can be written at X+ in tumbling co- 
ordinates Ha as, 



Wn 
Pn 
Jn 

{Ji)n 
{Ji)s 
Un 
Us 



Ws = Q, 
/3s = 0, 
Js = 0, 



3MA2ne-2«/M 



1-Cn 
■(1 + CjvCjv)' 

'{l + CsCs)' 



„/j,j 2i(^N cos KU+ (1 + Cat^) sin kuY 

(1 + CA^Civ)^ 

2„ /M 2«Cs cos KU + {1 + Cs^) sin 



(1 + CsCs)' 



(36) 



In this case we again set A/ = 1, A20 = 3 x 10^^ and 
K = 0.01. With this data, Jo shows 2'"^ order conver- 
gence to zero as for the Schwarzschild in tumbling co- 
ordinates. L2 norm of the error in the news also shows 
Qj.jjgj. convergence to zero as shown in figure (PJ and 
table ijirjl . The accuracy of news extraction remains very 
good as long as equator of one patch in Bondi coordinates 
doesn't go over the pole of the other patch in Bondi-Sachs 



convergence rate at u 


grids (a) & (b) 


(b) & (c) 


(c) & (d) 


4 


2.015 


2.002 


2.002 


8 


2.036 


1.982 


2.032 


12 


2.055 


1.943 


2.064 


16 


2.049 


1.906 


2.106 


20 


2.031 


1.862 


2.156 



TABLE II: Change of convergence rate with u for the L2 
norm of error in the news function for linearized Robinson- 
Trautman in tumbling coordinates for different grid resolu- 
tions in eg nation 1331 



coordinates due to tumbling. As such the accuracy of the 
news extraction does not depend on the initial amplitude 
of the waveform even it is increased or decreased by more 
that 10 orders of magnitude. 



VI. CONCLUSIONS 

We have developed a new scheme for extracting the 
news from characteristic numerical simulations of the 
spacetimes. Like our earlier scheme this scheme is based 
on a coordinate transformation from general Bondi-Sachs 
coordinates to Bondi coordinates. But in the present 
scheme the transformation is done in a different way. It 
is done so that we can easily use a uniform and constant 
Bondi angular grid for the news calculation. This has an 
advantage that while transforming coordinates one has 
to evolve ordinary differential equations instead of partial 
ones. Also the interpolation and angular differentiation 
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becomes easier and more accurate. We also introduce 
dU and dU as auxiliary variables in the code to over- 
come problems related to numerical differentiations. 

Like earlier schemes the present scheme also needs J — 
at I"*" on the initial time {u — 0) slice. It also needs 
the Bondi-Sachs and Bondi coordinates to match on the 
initial slice. We can also check the accuracy of coordinate 
transformation by monitoring values of Jq, which should 
be zero analytically. 

Our new scheme shows 2"'' order convergence of Jq to 
zero and overcomes the convergence problem of Jq faced 
by our earlier scheme for the Schwarzschild solution in 
tumbling coordinates. Also with this scheme the news 
for Schwarzschild solution in rotating/tumbling coordi- 
nates stays exactly zero, while with the earlier scheme it 
is not so small value which was converging to zero to 2"'^ 
order in grid size. In other cases also our new scheme 
shows 2"*^ order convergence in extracting the news. For 
the Robinson- Trautman test, it successfully extracts out 
the news even when the amplitude of the gravitational 
waves is very small and we are in a coordinate system 
which is rotating/tumbling with respect to the Bondi co- 
ordinates. The news calculation remains accurate even 
though the signal goes down many orders of magnitude 
with time. The relative accuracy of the news extraction 



does not depend on the initial amplitude of the waveform 
even when it is increased or decreased by more that 10 
orders of magnitude from the cases studied in this pa- 
per. This indicates our scheme gives very accurate and 
faithful results and the scheme is quite robust. 

Introducing dU and (5U as auxiliary variables seems 
to play a crucial role in getting the news accurately. In 
most cases, whichever scheme we may choose to calculate 
the news, it will indirectly contain at least third angular 
derivatives of U. Higher angular derivatives are known 
to create problems in the null code. To take care of this 
problem earlier three complex variables corresponding to 
3/3, 9 J and dK were introduced in j^. In a similar 
manner we suggest to introduce eth derivatives of U as 
new variables. It will be particularly useful when one is 
extracting gravitational waves from the simulations. All 
the tests we have done show significant improvement if 
we introduce these new variables. 

The accuracy of our new scheme is very encouraging 
and we expect that the present method will be useful 
for extracting gravitational radiation from more realistic 
physical systems/simulations. It could play an impor- 
tant role in Cauchy-Characteristic extraction. We plan 
to explore this direction in a future work. 



APPENDIX A: EQUATION FOR Ji 

The expression can be written as, 

= {-4duKou{{dX{Jo - Ko)P + {~Jq + Ko)P{dX + PBxX) + 2iPBuJm{UoiJQ ~ Ko)) 
+PduKoUo){-2P{Ui -Ui + AuLoiUu - Uu)) + 2iXm{ALjP{-W + W + PUaX))) 
+ (-P{Jo + Ko){dX + PBxX) - PgX(Jo + Kq) + 2Pdune{Uo{Jo + Kq))) 
{{2P{{Ui + Ue) + Au^{Uu + Uu))) + 2ne{AwP{W + W + PUoX)))) 
-{8uj{{dX + BxPX)^KqP^ + P^{^u)^Uo{KoUq + JqUq) + {^X + BxPX)P^QXJq 
-PPdu{dX + BxPX){2Ko(Jo + JaUa) + PduP{-dX JoUo)){P{Je + A„ Jo,n^) 

+ ^{A{2dJ + (1 - i)JoPX - (1 + i)JoPX) + A{2dJ + {1 + i)JoPX - (1 - i)JoPx))))/P 
-{8uj({dX)^KoP^ + {dX + BxPX)JoP{5XP - PduUa) + PduUo{Pdu{JaUo + KqUq)) 

+dX{PduPi~JoUa - 2KoUo)))iP{Ji + A„ Jo,„ lu) + '^{A{2dJ ~ (1 + i)JoPX + (1 - i)JoPx) 

+A{2dJ~ (1 - i)JoPX + {1 + i)JoPX))))/P + 16KQ{P^du{2e^'^du - dAuUj^{JoU^ + Uo{2KoUo + JqUq))) 

-Lo^P^ddA + ABxP)(pXJo + BxJqPX + dXKo) + dA{dXjQ + dXKo + BxKqPx)) 

+Pu;^PiiduiQA + BxPA){JoUo + KqUo) + dA{KoUo + JoU„)) + 9A„((gX + BxPx)iJ„Uo + K„U„) 

+QXiKoUo + JoUo))) + P^e^>'iBu)^ujVai))/il6P^Ko) (Al) 

I 
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